home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
The Datafile PD-CD 1 Issue 2
/
PDCD-1 - Issue 02.iso
/
_utilities
/
utilities
/
001
/
meschach
/
!Meschach
/
c
/
init
< prev
next >
Wrap
Text File
|
1994-01-13
|
6KB
|
299 lines
/**************************************************************************
**
** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
**
** Meschach Library
**
** This Meschach Library is provided "as is" without any express
** or implied warranty of any kind with respect to this software.
** In particular the authors shall not be liable for any direct,
** indirect, special, incidental or consequential damages arising
** in any way from use of the software.
**
** Everyone is granted permission to copy, modify and redistribute this
** Meschach Library, provided:
** 1. All copies contain this copyright notice.
** 2. All modified copies shall carry a notice stating who
** made the last modification and the date of such modification.
** 3. No charge is made for this software or works derived from it.
** This clause shall not be construed as constraining other software
** distributed on the same medium as this software, nor is a
** distribution fee considered a charge.
**
***************************************************************************/
/*
This is a file of routines for zero-ing, and initialising
vectors, matrices and permutations.
This is to be included in the matrix.a library
*/
static char rcsid[] = "$Id: init.c,v 1.6 1994/01/13 05:36:58 des Exp $";
#include <stdio.h>
#include "matrix.h"
/* v_zero -- zero the vector x */
VEC *v_zero(x)
VEC *x;
{
if ( x == VNULL )
error(E_NULL,"v_zero");
__zero__(x->ve,x->dim);
/* for ( i = 0; i < x->dim; i++ )
x->ve[i] = 0.0; */
return x;
}
/* iv_zero -- zero the vector ix */
IVEC *iv_zero(ix)
IVEC *ix;
{
int i;
if ( ix == IVNULL )
error(E_NULL,"iv_zero");
for ( i = 0; i < ix->dim; i++ )
ix->ive[i] = 0;
return ix;
}
/* m_zero -- zero the matrix A */
MAT *m_zero(A)
MAT *A;
{
int i, A_m, A_n;
Real **A_me;
if ( A == MNULL )
error(E_NULL,"m_zero");
A_m = A->m; A_n = A->n; A_me = A->me;
for ( i = 0; i < A_m; i++ )
__zero__(A_me[i],A_n);
/* for ( j = 0; j < A_n; j++ )
A_me[i][j] = 0.0; */
return A;
}
/* mat_id -- set A to being closest to identity matrix as possible
-- i.e. A[i][j] == 1 if i == j and 0 otherwise */
MAT *m_ident(A)
MAT *A;
{
int i, size;
if ( A == MNULL )
error(E_NULL,"m_ident");
m_zero(A);
size = min(A->m,A->n);
for ( i = 0; i < size; i++ )
A->me[i][i] = 1.0;
return A;
}
/* px_ident -- set px to identity permutation */
PERM *px_ident(px)
PERM *px;
{
int i, px_size;
u_int *px_pe;
if ( px == PNULL )
error(E_NULL,"px_ident");
px_size = px->size; px_pe = px->pe;
for ( i = 0; i < px_size; i++ )
px_pe[i] = i;
return px;
}
/* Pseudo random number generator data structures */
/* Knuth's lagged Fibonacci-based generator: See "Seminumerical Algorithms:
The Art of Computer Programming" sections 3.2-3.3 */
#ifdef ANSI_C
#ifndef LONG_MAX
#include <limits.h>
#endif
#endif
#ifdef LONG_MAX
#define MODULUS LONG_MAX
#else
#define MODULUS 1000000000L /* assuming long's at least 32 bits long */
#endif
#define MZ 0L
static long mrand_list[56];
static int started = FALSE;
static int inext = 0, inextp = 31;
/* mrand -- pseudo-random number generator */
#ifdef ANSI_C
double mrand(void)
#else
double mrand()
#endif
{
long lval;
static Real factor = 1.0/((Real)MODULUS);
if ( ! started )
smrand(3127);
inext = (inext >= 54) ? 0 : inext+1;
inextp = (inextp >= 54) ? 0 : inextp+1;
lval = mrand_list[inext]-mrand_list[inextp];
if ( lval < 0L )
lval += MODULUS;
mrand_list[inext] = lval;
return (double)lval*factor;
}
/* mrandlist -- fills the array a[] with len random numbers */
void mrandlist(a, len)
Real a[];
int len;
{
int i;
long lval;
static Real factor = 1.0/((Real)MODULUS);
if ( ! started )
smrand(3127);
for ( i = 0; i < len; i++ )
{
inext = (inext >= 54) ? 0 : inext+1;
inextp = (inextp >= 54) ? 0 : inextp+1;
lval = mrand_list[inext]-mrand_list[inextp];
if ( lval < 0L )
lval += MODULUS;
mrand_list[inext] = lval;
a[i] = (Real)lval*factor;
}
}
/* smrand -- set seed for mrand() */
void smrand(seed)
int seed;
{
int i;
mrand_list[0] = (123413*seed) % MODULUS;
for ( i = 1; i < 55; i++ )
mrand_list[i] = (123413*mrand_list[i-1]) % MODULUS;
started = TRUE;
/* run mrand() through the list sufficient times to
thoroughly randomise the array */
for ( i = 0; i < 55*55; i++ )
mrand();
}
#undef MODULUS
#undef MZ
#undef FAC
/* v_rand -- initialises x to be a random vector, components
independently & uniformly ditributed between 0 and 1 */
VEC *v_rand(x)
VEC *x;
{
/* int i; */
if ( ! x )
error(E_NULL,"v_rand");
/* for ( i = 0; i < x->dim; i++ ) */
/* x->ve[i] = rand()/((Real)MAX_RAND); */
/* x->ve[i] = mrand(); */
mrandlist(x->ve,x->dim);
return x;
}
/* m_rand -- initialises A to be a random vector, components
independently & uniformly distributed between 0 and 1 */
MAT *m_rand(A)
MAT *A;
{
int i /* , j */;
if ( ! A )
error(E_NULL,"m_rand");
for ( i = 0; i < A->m; i++ )
/* for ( j = 0; j < A->n; j++ ) */
/* A->me[i][j] = rand()/((Real)MAX_RAND); */
/* A->me[i][j] = mrand(); */
mrandlist(A->me[i],A->n);
return A;
}
/* v_ones -- fills x with one's */
VEC *v_ones(x)
VEC *x;
{
int i;
if ( ! x )
error(E_NULL,"v_ones");
for ( i = 0; i < x->dim; i++ )
x->ve[i] = 1.0;
return x;
}
/* m_ones -- fills matrix with one's */
MAT *m_ones(A)
MAT *A;
{
int i, j;
if ( ! A )
error(E_NULL,"m_ones");
for ( i = 0; i < A->m; i++ )
for ( j = 0; j < A->n; j++ )
A->me[i][j] = 1.0;
return A;
}
/* v_count -- initialises x so that x->ve[i] == i */
VEC *v_count(x)
VEC *x;
{
int i;
if ( ! x )
error(E_NULL,"v_count");
for ( i = 0; i < x->dim; i++ )
x->ve[i] = (Real)i;
return x;
}